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Abstract 

We demonstrate that the spatial resolution of images in optical tomography is not limited to 
the fundamental length scale of one transport mean free path. This result is facilitated by the 
introduction of novel corrections to the standard integral equations of scattering theory within the 
diffusion approximation to the radiative transport equation. 



1 



There has been considerable recent interest in the development of optical methods for 
tomographic imaging Q|. The physical problem that is considered is to recover the optical 
properties of the interior of an inhomogeneous medium from measurements taken on its 
surface. The starting point for the mathematical formulation of this inverse scattering 
problem (ISP) is a model for the propagation of light, typically taken to be the diffusion 
approximation (DA) to the radiative transport equation (RTE). The DA is valid when the 
energy density of the optical field varies slowly on the scale of the transport mean free path 
C.* . The DA breaks down in optically thin layers, near boundary surfaces, or near the source. 
One or more of these conditions are encountered in biomedical applications such as imaging 
of small animals 0] or of functional activity in the brain. 



Within the accuracy o: 



and analytic solutions 



;he DA, reconstruction algorithms based on both numerical 



5|, la] to the ISP have been described. Regardless of the method 
of inversion, the spatial resolution of reconstructed images is expected to be limited to 
t. This expectation is due to the intertwined effects of the ill-posedness of the ISP 
and intrinsic inaccuracies of the DA [3]. In this Letter, we introduce novel corrections to 
the integral equations of scattering theory within the DA. Using this result, we report the 
reconstruction of superresolved images whose spatial resolution is less than i*. 

We begin by considering the propagation of multiply-scattered light in an inhomogenous 
medium characterized by an absorption coefficient ;Ua(r). In what follows, we will neglect the 
contribution of ballistic photons and consider only diffuse photons whose specific intensity 
J(r, s) at the point r in the direction s is taken to obey the time- independent RTE 

s • V/(r, s) + (/i, + /i,)/(r, s) - /i, 1 ^^^'^(s, s')/(r, s') = ^(r, s) , (1) 

where fig is the scattering coefficient, A{s, s') is the scattering kernel, and S{r, s) is the 
source. The change in specific intensity due to spatial fiuctuations in /ia(r) can be obtained 
from the integral equation 

0(ri, Si; ra, §2) = J d^rd^ sG{yi, §1; r, s)G(r, s; rg, S2)5/ia(r) . (2) 

Here the data function 0(ri, §1; r2, §2) is proportional, to lowest order in ^/i^, to the change 
in specific intensity relative to a reference medium with absorption G is the Green's 
function for ((H) with = fJ'a^ ^/^a(i') = Aia(i') ~ /^a; I"!, Si and r2, §2 denote the position and 
direction of a unidirectional point source and detector, respectively. 



We now show that the integral equation Q may be used to obtain corrections to the 
usual formulation of scattering theory within the DA. To proceed, we note that, following 
Ref. the Green's function G{r, s; r', s') may be expanded in angular harmonics of s and 



s': 



G(r, s; r', s') = ^ (1 + Ts ■ Vr) (1 - ts' ■ VrO G(r, r') , (3) 

where i* = l/[fia + f^'s] with = (1 — g)fis, 9 being the anisotropy of the scattering kernel 
A. The Green's function G{r, r') satisfies the diffusion equation (— Z^qV^ + ao) G{r, r') = 
6{r — r'), where the diffusion coefficient Dq = l/3ci* and ao = c/i°. In addition, the Green's 
function must satisfy boundary conditions on the surface of the medium (or at infinity in 
the case of free boundaries). In general we will consider boundary conditions of the form 
G{r, r') +£n- VG(r, r') = 0, where fi is the outward unit normal to the surface bounding the 
medium and i is the extrapolation distance. Making use of (jHI) and performing the angular 
integration over s in Q we obtain 



(ri,Si;r2,S2) = ^^1^2 J d^r 



■*2 

'fc(r) , (4) 



^(ri, r)G'(r, ra) - — VrG'(ri, r) ■ VrG(r, rs 



where the differential operators = 1 — (— l)^£*Sfc ■ Vr^. with k = 1,2 and 6a = cSfia- Note 
that if the source and detector are oriented in the inward and outward normal directions, 
respectively, then (jH) becomes 



ri,-n(ri);r2,n(r2)) = ^ ^1 + J d^r 



G(ri,r)G(r,r2) 



-— VrG'(ri,r)- VrG'(r,r2) 



5a(r) , (5) 



where we have used the boundary conditions on G to evaluate the action of the operators. 
Eq. is the main theoretical result of this Letter. It may be viewed as providing corrections 
to the DA since the first term on the right hand side of (jSj) corresponds to the standard DA 
in an inhomogeneous absorbing medium. We note that the second term may be interpreted 
as defining an effective diffusion coefficient D{r) = Dq — (£*^/3)(5a(r) since the expression 
VrG(ri,r) ■ VrG'(r, r2) defines the diffusion kernel in a medium with an inhomogeneous 
diffusion coefficient jsj. 

For the remainder of this paper we will work in the planar measurement geometry, often 
encountered in small-animal imaging. In this case, becomes 

(f){hoi,rho2) = J d^rK{pi,p2;r)Sa{r) , (6) 



where pi denotes the transverse coordinates of a point source in the plane z = 0, p2 denotes 
the transverse coordinates of a point detector in the plane z = L, and the dependence of (p 
on Si and §2 is not made explicit. Evidently, from considerations of invariance of the kernel 
K{pi, P2',r) under translations of its transverse arguments, it can be seen that K may be 
expressed as the Fourier integral 

K{pi, P2] r) = -^^^ j d^qid^q2n{(li, q2, z) exp [z(qi - q2) ■ p - «(qi ■ Pi - q2 ■ P2)] , 

where r = {p,z). The function k may be obtained from the plane- wave expansion of the 
diffusion Green's function obeying appropriate boundary conditions. In the case of free 
boundaries, it is readily seen that k, is given by the expression 



167rD2g(qi)g(q2) 



fl*2 



1 + Y (Q(qi)Q(q2) - qi ■ q2) 



X 



l + r(Q(qi)+Q(q2))+r2g(qi)g(q2) exp [-g(qi)k| - g(q2)k - ^|] , (7) 



where Q(q) = (g^ + ao/DoY^"^ and we have assumed that §1 = §2 = z. 

Inversion of the integral equation ® may be carried out by analytic methods. These 
methods have been shown to be computationally efficient and may be applied to data sets 
consisting of a very large number of measurements j^O]. The approach taken is to construct 
the singular value decomposition of the linear operator K in the proper Hilbert space setting 
and then use this result to obtain the pseudoinverse solution to In this manner, it is 
possible to account for the effects of sampling and thereby obtain the best (in the sense of 
minimizing the appropriate norm) bandlimited approximation to 6a. Here we use this 
approach to simulate the reconstruction of a point absorber located at a point ro between 
the measurement planes with 6a{r) = A6{r — fq) for constant A. In this situation it is 
possible to calculate the data function within radiative transport theory, thus avoiding 
"inverse crime." To proceed, we require the Green's function G{r, s; r', s') for the RTE in 
a homogeneous infinite medium which we obtain as described in Ref. j^. Note that in this 
case, the angular integration over s in Q may be carried out analytically. 

The effects of corrections to the DA were studied in numerical simulations following the 
methods of Ref. Q]. The simulations were performed for a medium with optical properties 
similar to breast tissue in the near infrared j^. The background absorption and reduced 
scattering coefficients were given by /i° = 0.03 cm~^ and p'^ = 10 cm~^. The scattering 
kernel was taken to be of Henyey-Greenstein type with y4(s, s') = J^iLoQ^Pei^ ■ s') and 



g = 0.98. This choice of parameters corresponds to = 1 mm and Dq = 0.8 cm^ns"^. The 
separation between the measurement planes L was varied in order to explore the effects of 
the corrections. A single point absorber was placed at the midpoint of the measurement 
planes with fq = (0,0, L/2) and A = 1 cm^ ns^^. The sources and detectors were located 
on a square lattice with spacing h. The total number of source-detector pairs was varied, 
along with h, as indicated below. To demonstrate the stability of the reconstruction in the 
presence of noise, Gaussian noise of zero mean was added to the data at the 1% level, relative 
to the average signal. Note that the level of regularization was chosen to be the same for all 
reconstructions. 

Reconstruction of 6a{r) for a point absorber defines the point spread function (PSF) of the 
reconstruction algorithm. The resolution Ax is defined as the half width at half maximum 
of the PSF. In Fig. 1(a) we consider the case of a thick layer with L = 6.6 cm. The above 
parameters were chosen to be /i = 0.83 mm and A^ = 1.5 x 10^. PSFs with and without 
corrections are shown. It may be seen that the effect of the corrections is negligible in the 
case of a thick layer and the resolution Ax = 3.5i*. For the case of a layer of intermediate 
thickness with L = 1.1 cm, as shown in Fig. 1(b), the corrections have a more significant 
effect. In particular, with h = 0.28 mm and A^ = 1.2 X 10", we found that Ax = OM* 
for the uncorrected reconstruction and Ax = 0.7i* for the corrected reconstruction. The 
corrections are most significant for the thinnest layer considered in this study, achieving a 
factor of two improvement in resolution when L = 0.55 cm. In this case, with h = 0.14 mm 
and A^ = 1.9 X 10^^, we found that Ax = OAi* for the uncorrected case and Ax = 0.2^* 
for the corrected case as shown in Fig. 1(c). Note that the number of source-detector pairs 
A^ 10^ — 10^^ may be achieved in modern non-contact optical tomography systems |lO|. 

In conclusion, we have described a series of corrections to the usual formulation of the 
DA in optical tomography. We have found that these corrections give rise to superresolved 
images with resolution below i*. Several comments on these results are necessary. First, 
the effects of corrections were demonstrated to be most significant in optically thin layers. 
However, corrections to the DA may also be expected to be important for thick layers when 
inhomogeneities in the absorption are located near the surface. Second, the results of this 
study were obtained without resorting to so-called inverse crime. That is, forward scattering 
data was obtained from the full RTE under conditions when the DA is known to break down. 
Third, the use of analytic reconstruction algorithms was essential for handling the extremely 
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large data employed in this study. Finally, we note that higher order corrections to the DA 
are also expected to be important for the nonlinear ISP. 
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Fig. 1. Reconstructions of a point absorber for different thicknesses of the slab using the 
corrected (sohd curve) and uncorrected (dashed curve) DA. 
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